水稻微生物组时间序列分析2b-散点图拟合

之前分享了3月底发表的的 《水稻微生物组时间序列分析》的文章，大家对其中图绘制过程比较感兴趣。一上午收到了超30条留言，累计收到41个小伙伴的留言求精讲。

我们将花时间把此文的原始代码整理并精讲，祝有需要的小伙伴能有收获。

本系列按原文4幅组图，共分为4节。本文是第二节b，散点图拟合。

时间序列分析中，使用拟合曲线可以非常好的展示时间序列上组内数值波动较大散点图中潜在的规律。本文采用散点图，分别对水稻一生群落结果变化进行分析，同时对不同地点间群落距离变化进行分析，提示水稻根系微生物组以及不同品种和地点间随时间变化的规律。

前情提要

先了解一下图2的内容。

图2. 微生物组随时间变化的规律

图2. 田间水稻根系微生物组在8-10周趋于稳定。

A-D. 对两个水稻品种分别在两地进行的连续微生物组调查结果相关分析，发现8-10周后群落结构趋于稳定。

E. 所有时间点距离水稻最后取样点的Bray-Curtis距离。发现土壤呈现小幅波动变化，而根系呈现出先快后慢，逐渐趋近的变化规律。

F. 不同水稻品种在两个地点间的距离变化，发现土壤差异稳定，而根系微生物组差异随时间增长而趋于一致。

方法简介：A-D采用R的cor()函数计算pearson相关系数，并使用Corrplot包展示，时间轴使用pheatmap绘制热图展示。

E-F基于vegan包计算的所有样品两两间Bray-Curtis距离。分别挑选距离终点的距离，和两地间的距离与时间序列上的关系，并采用ggplot2可视化散点图，并添加拟合曲线方便观察变化规律。

图2. E-F散点图拟合展示距离变化

图2. E-F两个图看似很简单，其实说起来绘制过程还是非常复杂的。

其中图E是观察水稻一生菌群结构的变化，参照点为最后一次取样119天。按土壤/品种X地点共6个大组分别计算每个时间点距离终点均值Bray-Cutis距离。

图F则是为了说明植物根系的菌群结构受宿主调控影响，随时间逐渐成熟过程中两个地点间的距离会越来越少，即宿主要菌群结构中控制力逐渐增加。因为计算了两个种植地点间各品种和土壤间的距离，并散点图展示及拟合。来说明发现的结果。

图E有6个组，分别筛选样品、求均值，计算距离和绘图的代码为例进行讲解。

清空工作环境和加载包

# Set working enviroment in Rstudio, select Session - Set working directory - To source file location, default is runing directory rm(list=ls()) # clean enviroment object library("ggplot2") # load related packages library("grid") library("scales") library("vegan") library("agricolae") library("ggbiplot") library("dplyr") library("ggrepel") # Set ggplot2 drawing parameter, such as axis line and text size, lengend and title size, and so on. main_theme = theme(panel.background=element_blank(),                    panel.grid=element_blank(),                    axis.line.x=element_line(size=.5, colour="black"),                    axis.line.y=element_line(size=.5, colour="black"),                    axis.ticks=element_line(color="black"),                    axis.text=element_text(color="black", size=7),                    legend.position="right",                    legend.background=element_blank(),                    legend.key=element_blank(),                    legend.text= element_text(size=7),                    text=element_text(family="sans", size=7))

读入实验设计和OTU表

# Public file 1. "design.txt"  Design of experiment design = read.table("../data/design.txt", header=T, row.names= 1, sep="\t") # setting subset design if (TRUE){    design = subset(design,groupID %in% c("A50Cp0","A50Cp1","A50Cp10","A50Cp112","A50Cp119","A50Cp14","A50Cp2","A50Cp21","A50Cp28","A50Cp3","A50Cp35","A50Cp42","A50Cp49","A50Cp63","A50Cp7","A50Cp70","A50Cp77","A50Cp84","A50Cp91","A50Cp98","A50Sz0","A50Sz1","A50Sz10","A50Sz118","A50Sz13","A50Sz2","A50Sz20","A50Sz27","A50Sz3","A50Sz34","A50Sz41","A50Sz48","A50Sz5","A50Sz56","A50Sz62","A50Sz69","A50Sz7","A50Sz76","A50Sz83","A50Sz90","A50Sz97","IR24Cp0","IR24Cp1","IR24Cp10","IR24Cp112","IR24Cp119","IR24Cp14","IR24Cp2","IR24Cp21","IR24Cp28","IR24Cp3","IR24Cp35","IR24Cp42","IR24Cp49","IR24Cp63","IR24Cp7","IR24Cp70","IR24Cp77","IR24Cp84","IR24Cp91","IR24Cp98","IR24Sz0","IR24Sz1","IR24Sz10","IR24Sz118","IR24Sz13","IR24Sz2","IR24Sz20","IR24Sz27","IR24Sz3","IR24Sz34","IR24Sz41","IR24Sz48","IR24Sz5","IR24Sz56","IR24Sz62","IR24Sz69","IR24Sz7","IR24Sz76","IR24Sz83","IR24Sz90","IR24Sz97","soilCp0","soilCp42","soilCp49","soilCp57","soilCp63","soilCp77","soilCp83","soilCp91","soilCp98","soilCp118","soilSz0","soilSz41","soilSz48","soilSz54","soilSz62","soilSz76","soilSz84","soilSz90","soilSz97","soilSz119") ) # select group1 } otu_table = read.delim("../data/otu_table.txt", row.names= 1,  header=T, sep="\t")

编写计算每种类型(指定范围)至距离终点的函数，方便简化代码和重用性

calc_dis = function(subset){ # 筛选名中包含部分关键字的组 idx = grepl(subset, design$groupID) sub_design = design[idx,] print(paste("Number of group: ",length(unique(sub_design$group)),sep="")) # show group numbers # 筛选样品，计算距离矩阵 # 有些样品异常会剔除，需要交叉筛选 idx=rownames(sub_design) %in% colnames(otu_table) sub_design=sub_design[idx,] count = otu_table[, rownames(sub_design)] norm = t(t(count)/colSums(count,na=T)) * 100 # normalization to total 100 norm=as.data.frame(norm) # 筛选119天，作为终点，再计算各点到终点的距离，观察群落结构变化 final=rownames(sub_design[sub_design$day2 %in% 119,]) ck=norm[,final] # 筛选并计算均值 ck_mean= rowMeans(ck) norm$ck_mean=ck_mean # 均值添加到OTU表 # 计算包括终点均值的所有样品bray距离 bray_curtis = vegdist(t(norm), method = "bray") bray_curtis= as.matrix(bray_curtis) # 计算各组内部差异程度 # 建立一个存储组内差异的数据框 dat = t(as.data.frame(c("sampleA","sampleB","0.15","group","genosite"))) colnames(dat) = c("sampleA","sampleB","distance","group","type") rownames(dat) = c("test") # 每个样品与final对应的距离 for (i in sort(unique(sub_design$day2))){  group = rownames(sub_design[sub_design$day2 %in% i,])  for (m in 1:(length(group)-1)) {    x = c(group[m],"ck_mean",bray_curtis[group[m],"ck_mean"],i,subset)    dat=rbind(dat,x)  } } dat = as.data.frame(dat[-1,], stringsAsFactors=F) # 删除首行框架数据 # dat = dat[dat$distance != 0,] # # 距离转换为3位数值 dat$distance=round(as.numeric(as.character(dat$distance)), digits=3) # 分组添加levels，分组有时间顺序 dat$group = factor(dat$group, levels=unique(dat$group)) return(dat) }

计算第一大类A50品种在Cp地点下每个点到终点的距离

dat = calc_dis("A50Cp")

计算第一次结果赋给all，以后每类别修改品种和地点重新运行上面55 - 103行后，cbind追加dat至all中

all = dat

分别按”A50”,”IR24”, X “Cp”,”Sz”筛选数据

for (i in c("A50Sz", "IR24Cp", "IR24Sz", "soilCp", "soilSz")){  dat = calc_dis(i)  all = rbind(all, dat) }

先用各组箱线图查看数据分布

p = ggplot(all, aes(x=group, y=distance, color=type)) +  geom_boxplot(alpha=1, outlier.size=0, size=0.7, width=0.5, fill="transparent") +  labs(x="Groups", y="Bray-Curtis distance") + theme_classic()   p=p+geom_jitter( position=position_jitter(0.17), size=1, alpha=0.7) p=p+theme(axis.text=element_text(angle=45,vjust=1, hjust=1)) p ggsave(paste("e.boxplot",".pdf", sep=""), p, width = 8, height =5 )

我们已经能看到各品种距离下降，而土壤中变化不大的趋势。使用散点图+拟合曲线应该可以更直观。

开始尝试散点图拟合

require(splines) require(MASS) # 为什么拟合曲线需要X轴为数字，因子转换为数字需要先转字符再转换数值，否则会无法生成拟合曲线 all$group = as.numeric(as.character(all$group)) p = ggplot(all, aes(x=group, y=distance, color=type)) +  geom_point(alpha=1, size=0.7) + geom_jitter()+  labs(x="Groups", y="Bray-Curtis distance") +  theme(axis.text=element_text(angle=45,vjust=1, hjust=1)) +  # 拟合只有loess方法初始为曲线，lm方法默认为直线，但可以添加公式  main_theme + geom_smooth(method = "lm", formula = y ~ poly(x,3)) # , formula = y ~ ns(x,3) p ggsave(paste("e.beta_scatter_bray_curtis_mean6",".pdf", sep=""), p, width = 4, height =2.5 )

我们的另一个备选方案，拆线图也很好看

library(ggpubr) p=ggline(all, x="group", y="distance", add = "mean_se",  color = "type",       palette = "jco",legend = "right") +main_theme p ggsave(paste("e.line",".pdf", sep=""), p, width = 5, height = 3)

对于图F的方法与E类似，仅是从任意时间点到终点的距离，变为了同一品种同一时间两地点间距离。具体代码见下面链接，就不在这里一一赘述了。

本分析的全部文件和代码，会在 https://github.com/YongxinLiu/RiceTimeCourse 上持续更新，也可以后台回复“时间序列”获得百度云下载链接。

如果本文分享的技术帮助了你的科研，欢迎引用下文，支持国产SCI越来越好。

Citition: Zhang, J., Zhang, N., Liu, Y.X., Zhang, X., Hu, B., Qin, Y., Xu, H., Wang, H., Guo, X., Qian, J., et al. (2018). Root microbiota shift in rice correlates with resident time in the field and developmental stage. Sci China Life Sci 61, https://doi.org/10.1007/s11427-018-9284-4

猜你喜欢

系列教程：微生物组入门 Biostar 微生物组  宏基因组

专业技能：生信宝典 学术图表 高分文章 不可或缺的人

一文读懂：宏基因组 寄生虫益处 进化树

必备技能：提问 搜索  Endnote

文献阅读 热心肠 SemanticScholar Geenmedical

扩增子分析：图表解读 分析流程 统计绘图

在线工具：16S预测培养基 生信绘图

科研经验：云笔记  云协作 公众号

编程模板 Shell  R Perl

生物科普  生命大跃进  细胞暗战 人体奥秘  

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外1700+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决群内讨论，问题不私聊，帮助同行。

学习16S扩增子、宏基因组科研思路和分析实战，关注“宏基因组”

点击阅读原文，跳转最新文章目录阅读